# Replication File for
# "Behavioral Consequences of Open Candidate Recruitment"
# Authors: Jochen Rehmert


### set directory, read data, load packages
setwd("...")

# libraries
library(foreign)
library(ggplot2);library(scales)
library(pscl);library(countreg);library(texreg)
library(sandwich);library(lmtest);library(stargazer);library(MASS)

# generich functions
char = function(x){as.character(x)}
num = function(x){as.numeric(char(x))}

kobo.parties <- c("DPJ","LDP","YP","JRP","NPD")

dat <- read.dta("kobo_activity.dta")
dat <- dat[dat$speaker == 0 & dat$replacement == 0 & dat$drop_out == 0,]

dat$party_id <- as.character(dat$party_id)
dat$LDP <- as.numeric(dat$party_id == "LDP")
dat$DPJ <- as.numeric(dat$party_id == "DPJ")
dat$SDP <- as.numeric(dat$party_id == "SDP")
dat$YP <- as.numeric(dat$party_id == "YP")
dat$JRP <- as.numeric(dat$party_id == "JRP")

dat$margin <- as.numeric(as.character(gsub(",",".",dat$margin)))

dat$kobo_f <- NA
dat$kobo_f[dat$kobo == 1 & dat$term <= 1] <- 1
dat$kobo_f[dat$kobo != 1 ] <- 0
dat$kobo_f[dat$kobo == 1 & dat$term >= 2] <- 0


###############################################################################
#                               Table A2                                      #
###############################################################################

summary(quest.negbin.b <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                   parliamentary_office + executive_office + margin + kobo  + government, 
                                 data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(quest.negbin.b$model)
quest.negbin.b.se <- coeftest(quest.negbin.b, vcov = vcovCL(quest.negbin.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(quest.negbin.sen.b <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                       parliamentary_office + executive_office + margin + kobo  * term + government, 
                                     data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(quest.negbin.sen.b$model)
quest.negbin.sen.b.se <- coeftest(quest.negbin.sen.b, vcov = vcovCL(quest.negbin.sen.b, cluster = dat[rows,"candidate_id"]))[,2]


summary(quest.negbin.gov.b <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                       parliamentary_office + executive_office + margin + kobo * government, 
                                     data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(quest.negbin.gov.b$model)
quest.negbin.gov.b.se <- coeftest(quest.negbin.gov.b, vcov = vcovCL(quest.negbin.gov.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(quest.negbin.r <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                   parliamentary_office + executive_office + margin + kobo_district  + government, 
                                 data = dat[which(dat$party_id %in% kobo.parties & dat$kobo_f == 0),]))
rows <- row.names(quest.negbin.r$model)
quest.negbin.r.se <- coeftest(quest.negbin.r, vcov = vcovCL(quest.negbin.r, cluster = dat[rows,"candidate_id"]))[,2]

summary(quest.negbin.gov.r <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                       parliamentary_office + executive_office + margin + kobo_district  * government, 
                                     data = dat[which(dat$party_id %in% kobo.parties & dat$kobo_f == 0),]))
rows <- row.names(quest.negbin.gov.r$model)
quest.negbin.gov.r.se <- coeftest(quest.negbin.gov.r, vcov = vcovCL(quest.negbin.gov.r, cluster = dat[rows,"candidate_id"]))[,2]

summary(quest.negbin.sen.r <- glm.nb(questions ~  as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                       parliamentary_office + executive_office + margin + kobo_district * term + government, 
                                     data = dat[which(dat$party_id %in% kobo.parties & dat$kobo_f == 0),]))
rows <- row.names(quest.negbin.sen.r$model)
quest.negbin.sen.r.se <- coeftest(quest.negbin.sen.r, vcov = vcovCL(quest.negbin.sen.r, cluster = dat[rows,"candidate_id"]))[,2]


stargazer(list(quest.negbin.b, quest.negbin.sen.b, quest.negbin.gov.b,
               quest.negbin.r, quest.negbin.sen.r, quest.negbin.gov.r),
          se = list(quest.negbin.b.se, quest.negbin.sen.b.se, quest.negbin.gov.b.se,
                    quest.negbin.r.se, quest.negbin.sen.r.se, quest.negbin.gov.r.se),
          omit = c("party_id","e2005","e2009","e2012"))

###############################################################################
#                               Table A3                                      #
###############################################################################

# Operationalization A
summary(quest.negbin.a <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                   parliamentary_office + executive_office + margin + kobo_f  + government, 
                                 data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(quest.negbin.a$model)
quest.negbin.a.se <- coeftest(quest.negbin.a, vcov = vcovCL(quest.negbin.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(quest.negbin.gov.a <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                       parliamentary_office + executive_office + margin + kobo_f  * government, 
                                     data = dat[which(dat$party_id %in% kobo.parties),])) 
rows <- row.names(quest.negbin.gov.a$model)
quest.negbin.gov.a.se <- coeftest(quest.negbin.gov.a, vcov = vcovCL(quest.negbin.gov.a, cluster = dat[rows,"candidate_id"]))[,2]

stargazer(list(quest.negbin.a,  quest.negbin.gov.a),
          se = list(quest.negbin.a.se, quest.negbin.gov.a.se))


###############################################################################
#                               Table A4                                      #
###############################################################################

summary(co.negbin.b <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                parliamentary_office + executive_office + margin + kobo  + government, 
                              data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(co.negbin.b$model)
co.negbin.b.se <- coeftest(co.negbin.b, vcov = vcovCL(co.negbin.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(co.negbin.sen.b <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                    parliamentary_office + executive_office + margin + kobo  * term + government, 
                                  data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(co.negbin.sen.b$model)
co.negbin.sen.b.se <- coeftest(co.negbin.sen.b, vcov = vcovCL(co.negbin.sen.b, cluster = dat[rows,"candidate_id"]))[,2]


summary(co.negbin.gov.b <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                    parliamentary_office + executive_office + margin + kobo * government, 
                                  data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(co.negbin.gov.b$model)
co.negbin.gov.b.se <- coeftest(co.negbin.gov.b, vcov = vcovCL(co.negbin.gov.b, cluster = dat[rows,"candidate_id"]))[,2]


summary(co.negbin.r <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                parliamentary_office + executive_office + margin + kobo_district  + government, 
                              data = dat[which(dat$party_id %in% kobo.parties & dat$kobo_f == 0),]))
rows <- row.names(co.negbin.r$model)
co.negbin.r.se <- coeftest(co.negbin.r, vcov = vcovCL(co.negbin.r, cluster = dat[rows,"candidate_id"]))[,2]


summary(co.negbin.gov.r <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                    parliamentary_office + executive_office + margin + kobo_district  * government, 
                                  data = dat[which(dat$party_id %in% kobo.parties & dat$kobo_f == 0),]))
rows <- row.names(co.negbin.gov.r$model)
co.negbin.gov.r.se <- coeftest(co.negbin.gov.r, vcov = vcovCL(co.negbin.gov.r, cluster = dat[rows,"candidate_id"]))[,2]

summary(co.negbin.sen.r <- glm.nb(num_pmb_cosponsored ~  as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                    parliamentary_office + executive_office + margin + kobo_district * term + government, 
                                  data = dat[which(dat$party_id %in% kobo.parties & dat$kobo_f == 0),]))
rows <- row.names(co.negbin.sen.r$model)
co.negbin.sen.r.se <- coeftest(co.negbin.sen.r, vcov = vcovCL(co.negbin.sen.r, cluster = dat[rows,"candidate_id"]))[,2]


stargazer(list(co.negbin.b, co.negbin.sen.b, co.negbin.gov.b,
               co.negbin.r, co.negbin.sen.r, co.negbin.gov.r),
          se = list(co.negbin.b.se, co.negbin.sen.b.se, co.negbin.gov.b.se,
                    co.negbin.r.se, co.negbin.sen.r.se, co.negbin.gov.r.se))

###############################################################################
#                               Table A5                                      #
###############################################################################

summary(co.negbin.a <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                parliamentary_office + executive_office + margin + kobo_f  + government, 
                              data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(co.negbin.a$model)
co.negbin.a.se <- coeftest(co.negbin.a, vcov = vcovCL(co.negbin.a, cluster = dat[rows,"candidate_id"]))[,2]

summary(co.negbin.gov.a <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                    parliamentary_office + executive_office + margin + kobo_f  * government, 
                                  data = dat[which(dat$party_id %in% kobo.parties),])) 
rows <- row.names(co.negbin.gov.a$model)
co.negbin.gov.a.se <- coeftest(co.negbin.gov.a, vcov = vcovCL(co.negbin.gov.a, cluster = dat[rows,"candidate_id"]))[,2]

stargazer(list(co.negbin.a,  co.negbin.gov.a),
          se = list(co.negbin.a.se, co.negbin.gov.a.se))

###############################################################################
#                               Table A6                                      #
###############################################################################
summary(ini.negbin.b <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                 parliamentary_office + executive_office + margin + kobo  + government, 
                               data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(ini.negbin.b$model)
ini.negbin.b.se <- coeftest(ini.negbin.b, vcov = vcovCL(ini.negbin.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(ini.negbin.sen.b <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + kobo  * term + government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(ini.negbin.sen.b$model)
ini.negbin.sen.b.se <- coeftest(ini.negbin.sen.b, vcov = vcovCL(ini.negbin.sen.b, cluster = dat[rows,"candidate_id"]))[,2]


summary(ini.negbin.gov.b <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + kobo * government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(ini.negbin.gov.b$model)
ini.negbin.gov.b.se <- coeftest(ini.negbin.gov.b, vcov = vcovCL(ini.negbin.gov.b, cluster = dat[rows,"candidate_id"]))[,2]


summary(ini.negbin.r <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                 parliamentary_office + executive_office + margin + kobo_district  + government, 
                               data = dat[which(dat$party_id %in% kobo.parties & dat$kobo_f == 0),]))
rows <- row.names(ini.negbin.r$model)
ini.negbin.r.se <- coeftest(ini.negbin.r, vcov = vcovCL(ini.negbin.r, cluster = dat[rows,"candidate_id"]))[,2]

summary(ini.negbin.sen.r <- glm.nb(num_pmbs_initiated ~  as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + kobo_district * term + government, 
                                   data = dat[which(dat$party_id %in% kobo.parties & dat$kobo_f == 0),]))
rows <- row.names(ini.negbin.sen.r$model)
ini.negbin.sen.r.se <- coeftest(ini.negbin.sen.r, vcov = vcovCL(ini.negbin.sen.r, cluster = dat[rows,"candidate_id"]))[,2]

summary(ini.negbin.gov.r <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + kobo_district  * government, 
                                   data = dat[which(dat$party_id %in% kobo.parties & dat$kobo_f == 0),]))
rows <- row.names(ini.negbin.gov.r$model)
ini.negbin.gov.r.se <- coeftest(ini.negbin.gov.r, vcov = vcovCL(ini.negbin.gov.r, cluster = dat[rows,"candidate_id"]))[,2]


stargazer(list(ini.negbin.b, ini.negbin.sen.b, ini.negbin.gov.b,
               ini.negbin.r, ini.negbin.sen.r, ini.negbin.gov.r),
          se = list(ini.negbin.b.se, ini.negbin.sen.b.se, ini.negbin.gov.b.se,
                    ini.negbin.r.se, ini.negbin.sen.r.se, ini.negbin.gov.r.se))

###############################################################################
#                               Table A7                                      #
###############################################################################
summary(ini.negbin.a <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + kobo_f  + government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),])) 
rows <- row.names(ini.negbin.a$model)
ini.negbin.a.se <- coeftest(ini.negbin.a, vcov = vcovCL(ini.negbin.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(ini.negbin.gov.a <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + kobo_f  * government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),])) 
rows <- row.names(ini.negbin.gov.a$model)
ini.negbin.gov.a.se <- coeftest(ini.negbin.gov.a, vcov = vcovCL(ini.negbin.gov.a, cluster = dat[rows,"candidate_id"]))[,2]

stargazer(list(ini.negbin.a, ini.negbin.gov.a),
          se = list(ini.negbin.a.se, ini.negbin.gov.a.se))

###############################################################################
#                               Table A8                                      #
###############################################################################

summary(quest.zero.a <- zeroinfl(questions ~  term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                   parliamentary_office  + margin + kobo_f  + government | -1 + kobo_f +  executive_office, data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(quest.zero.a$model)
quest.zero.a.se <- coeftest(quest.zero.a, vcov = vcovCL(quest.zero.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(quest.zero.gov.a <- zeroinfl(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                       parliamentary_office  + margin + kobo_f * government | - 1 + kobo_f + executive_office, data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(quest.zero.gov.a$model)
quest.zero.gov.a.se <- coeftest(quest.zero.gov.a, vcov = vcovCL(quest.zero.gov.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(quest.zero.b <- zeroinfl(questions ~  term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                   parliamentary_office  + margin + kobo  + government | kobo + executive_office , data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(quest.zero.b$model)
quest.zero.b.se <- coeftest(quest.zero.b, vcov = vcovCL(quest.zero.b, cluster = dat[rows,"candidate_id"]))[,2]


summary(quest.zero.gov.b <- zeroinfl(questions ~  term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                       parliamentary_office  + margin + kobo * government | kobo + executive_office, data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(quest.zero.gov.b$model)
quest.zero.gov.b.se <- coeftest(quest.zero.gov.b, vcov = vcovCL(quest.zero.gov.b, cluster = dat[rows,"candidate_id"]))[,2]


summary(quest.zero.sen.b <- zeroinfl(questions ~  term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                       parliamentary_office  + margin + kobo  * term + government | kobo + executive_office, data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(quest.zero.sen.b$model)
quest.zero.sen.b.se <- coeftest(quest.zero.sen.b, vcov = vcovCL(quest.zero.sen.b, cluster = dat[rows,"candidate_id"]))[,2]

texreg(list(quest.zero.a, quest.zero.gov.a,quest.zero.b, quest.zero.gov.b,quest.zero.sen.b),
       se = list(quest.zero.a.se, quest.zero.gov.a.se,quest.zero.b.se, quest.zero.gov.b.se, quest.zero.sen.b.se), 
       stars = c(0.01,0.05, 0.1), digits = 3)


###############################################################################
#                               Table A9                                      #
###############################################################################

summary(co.zero.a <- zeroinfl(num_pmb_cosponsored ~  term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                parliamentary_office  + margin + kobo_f  + government |  kobo_f +  executive_office, data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(co.zero.a$model)
co.zero.a.se <- coeftest(co.zero.a, vcov = vcovCL(co.zero.a, cluster = dat[rows,"candidate_id"]))[,2]

summary(co.zero.gov.a <- zeroinfl(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                    parliamentary_office  + margin + kobo_f * government |  kobo_f + executive_office, data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(co.zero.gov.a$model)
co.zero.gov.a.se <- coeftest(co.zero.gov.a, vcov = vcovCL(co.zero.gov.a, cluster = dat[rows,"candidate_id"]))[,2]

summary(co.zero.b <- zeroinfl(num_pmb_cosponsored ~  term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                parliamentary_office  + margin + kobo  + government | kobo + executive_office , data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(co.zero.b$model)
co.zero.b.se <- coeftest(co.zero.b, vcov = vcovCL(co.zero.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(co.zero.gov.b <- zeroinfl(num_pmb_cosponsored ~  term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                    parliamentary_office  + margin + kobo * government | kobo + executive_office, data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(co.zero.gov.b$model)
co.zero.gov.b.se <- coeftest(co.zero.gov.b, vcov = vcovCL(co.zero.gov.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(co.zero.sen.b <- zeroinfl(num_pmb_cosponsored ~  term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                    parliamentary_office  + margin + kobo  * term + government | kobo + executive_office, data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(co.zero.sen.b$model)
co.zero.sen.b.se <- coeftest(co.zero.sen.b, vcov = vcovCL(co.zero.sen.b, cluster = dat[rows,"candidate_id"]))[,2]

library(texreg)
texreg(list(co.zero.a, co.zero.gov.a,co.zero.b, co.zero.gov.b,co.zero.sen.b),
       override.se = list(co.zero.a.se, co.zero.gov.a.se,co.zero.b.se, co.zero.gov.b.se, co.zero.sen.b.se), 
       stars = c(0.01,0.05, 0.1), digits = 3)

###############################################################################
#                               Table A10                                     #
###############################################################################

summary(ini.zero.a <- zeroinfl(num_pmbs_initiated ~  term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                 parliamentary_office  + margin + kobo_f  + government |  kobo_f +  executive_office, data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(ini.zero.a$model)
ini.zero.a.se <- coeftest(ini.zero.a, vcov = vcovCL(ini.zero.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(ini.zero.gov.a <- zeroinfl(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office  + margin + kobo_f * government |  kobo_f + executive_office, data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(ini.zero.gov.a$model)
ini.zero.gov.a.se <- coeftest(ini.zero.gov.a, vcov = vcovCL(ini.zero.gov.a, cluster = dat[rows,"candidate_id"]))[,2]

summary(ini.zero.b <- zeroinfl(num_pmbs_initiated ~  term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                 parliamentary_office  + margin + kobo  + government | kobo + executive_office , data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(ini.zero.b$model)
ini.zero.b.se <- coeftest(ini.zero.b, vcov = vcovCL(ini.zero.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(ini.zero.gov.b <- zeroinfl(num_pmbs_initiated ~  term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office  + margin + kobo * government | kobo + executive_office, data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(ini.zero.gov.b$model)
ini.zero.gov.b.se <- coeftest(ini.zero.gov.b, vcov = vcovCL(ini.zero.gov.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(ini.zero.sen.b <- zeroinfl(num_pmbs_initiated ~  term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office  + margin + kobo  * term + government | kobo + executive_office, data = dat[which(dat$party_id %in% kobo.parties),], dist = "negbin", link = "logit"))
rows <- row.names(ini.zero.sen.b$model)
ini.zero.sen.b.se <- coeftest(ini.zero.sen.b, vcov = vcovCL(ini.zero.sen.b, cluster = dat[rows,"candidate_id"]))[,2]

texreg(list(ini.zero.a, ini.zero.gov.a,ini.zero.b, ini.zero.gov.b,ini.zero.sen.b),
       se = list(ini.zero.a.se, ini.zero.gov.a.se,ini.zero.b.se, ini.zero.gov.b.se, ini.zero.sen.b.se), 
       stars = c(0.01,0.05, 0.1), digits = 3)

###############################################################################
#                               Table A11                                     #
###############################################################################

summary(quest.negbin.a <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                 parliamentary_office + executive_office + margin + pure_smd + kobo_f  + government, 
                               data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(quest.negbin.a$model)
quest.negbin.a.se <- coeftest(quest.negbin.a, vcov = vcovCL(quest.negbin.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(quest.negbin.gov.a <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + pure_smd + kobo_f  * government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),])) 
rows <- row.names(quest.negbin.gov.a$model)
quest.negbin.gov.a.se <- coeftest(quest.negbin.gov.a, vcov = vcovCL(quest.negbin.gov.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(quest.negbin.b <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                 parliamentary_office + executive_office + margin + pure_smd + kobo  + government, 
                               data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(quest.negbin.b$model)
quest.negbin.b.se <- coeftest(quest.negbin.b, vcov = vcovCL(quest.negbin.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(quest.negbin.sen.b <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + pure_smd + kobo  * term + government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(quest.negbin.sen.b$model)
quest.negbin.sen.b.se <- coeftest(quest.negbin.sen.b, vcov = vcovCL(quest.negbin.sen.b, cluster = dat[rows,"candidate_id"]))[,2]


summary(quest.negbin.gov.b <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + pure_smd + kobo * government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(quest.negbin.gov.b$model)
quest.negbin.gov.b.se <- coeftest(quest.negbin.gov.b, vcov = vcovCL(quest.negbin.gov.b, cluster = dat[rows,"candidate_id"]))[,2]


stargazer(list(quest.negbin.a, quest.negbin.gov.a, quest.negbin.b, quest.negbin.sen.b, quest.negbin.gov.b),
          se = list(quest.negbin.a.se, quest.negbin.gov.a.se, quest.negbin.b.se, quest.negbin.sen.b.se, quest.negbin.gov.b.se),
          omit = "party_id")



###############################################################################
#                               Table A12                                     #
###############################################################################

summary(ini.negbin.a <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                 parliamentary_office + executive_office + margin + pure_smd + kobo_f  + government, 
                               data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(ini.negbin.a$model)
ini.negbin.a.se <- coeftest(ini.negbin.a, vcov = vcovCL(ini.negbin.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(ini.negbin.gov.a <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + pure_smd + kobo_f  * government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),])) 
rows <- row.names(ini.negbin.gov.a$model)
ini.negbin.gov.a.se <- coeftest(ini.negbin.gov.a, vcov = vcovCL(ini.negbin.gov.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(ini.negbin.b <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                 parliamentary_office + executive_office + margin + pure_smd + kobo  + government, 
                               data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(ini.negbin.b$model)
ini.negbin.b.se <- coeftest(ini.negbin.b, vcov = vcovCL(ini.negbin.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(ini.negbin.sen.b <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + pure_smd + kobo  * term + government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(ini.negbin.sen.b$model)
ini.negbin.sen.b.se <- coeftest(ini.negbin.sen.b, vcov = vcovCL(ini.negbin.sen.b, cluster = dat[rows,"candidate_id"]))[,2]


summary(ini.negbin.gov.b <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + pure_smd + kobo * government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(ini.negbin.gov.b$model)
ini.negbin.gov.b.se <- coeftest(ini.negbin.gov.b, vcov = vcovCL(ini.negbin.gov.b, cluster = dat[rows,"candidate_id"]))[,2]


stargazer(list(ini.negbin.a, ini.negbin.gov.a, ini.negbin.b, ini.negbin.sen.b, ini.negbin.gov.b),
          se = list(ini.negbin.a.se, ini.negbin.gov.a.se, ini.negbin.b.se, ini.negbin.sen.b.se, ini.negbin.gov.b.se))


###############################################################################
#                               Table A13                                     #
###############################################################################

summary(co.negbin.a <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                 parliamentary_office + executive_office + margin + pure_smd + kobo_f  + government, 
                               data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(co.negbin.a$model)
co.negbin.a.se <- coeftest(co.negbin.a, vcov = vcovCL(co.negbin.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(co.negbin.gov.a <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + pure_smd + kobo_f  * government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),])) 
rows <- row.names(co.negbin.gov.a$model)
co.negbin.gov.a.se <- coeftest(co.negbin.gov.a, vcov = vcovCL(co.negbin.gov.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(co.negbin.b <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                 parliamentary_office + executive_office + margin + pure_smd + kobo  + government, 
                               data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(co.negbin.b$model)
co.negbin.b.se <- coeftest(co.negbin.b, vcov = vcovCL(co.negbin.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(co.negbin.sen.b <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + pure_smd + kobo  * term + government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(co.negbin.sen.b$model)
co.negbin.sen.b.se <- coeftest(co.negbin.sen.b, vcov = vcovCL(co.negbin.sen.b, cluster = dat[rows,"candidate_id"]))[,2]


summary(co.negbin.gov.b <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + margin + pure_smd + kobo * government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(co.negbin.gov.b$model)
co.negbin.gov.b.se <- coeftest(co.negbin.gov.b, vcov = vcovCL(co.negbin.gov.b, cluster = dat[rows,"candidate_id"]))[,2]


stargazer(list(co.negbin.a, co.negbin.gov.a, co.negbin.b, co.negbin.sen.b, co.negbin.gov.b),
          se = list(co.negbin.a.se, co.negbin.gov.a.se, co.negbin.b.se, co.negbin.sen.b.se, co.negbin.gov.b.se))




###############################################################################
#                               Table A14                                     #
###############################################################################

summary(quest.negbin.a <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                   parliamentary_office + executive_office + zombie_now + kobo_f  + government, 
                                 data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(quest.negbin.a$model)
quest.negbin.a.se <- coeftest(quest.negbin.a, vcov = vcovCL(quest.negbin.a, cluster = dat[rows,"candidate_id"]))[,2]

summary(quest.negbin.gov.a <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                       parliamentary_office + executive_office + zombie_now + kobo_f  * government, 
                                     data = dat[which(dat$party_id %in% kobo.parties),])) 
rows <- row.names(quest.negbin.gov.a$model)
quest.negbin.gov.a.se <- coeftest(quest.negbin.gov.a, vcov = vcovCL(quest.negbin.gov.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(quest.negbin.b <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                   parliamentary_office + executive_office + zombie_now + kobo  + government, 
                                 data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(quest.negbin.b$model)
quest.negbin.b.se <- coeftest(quest.negbin.b, vcov = vcovCL(quest.negbin.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(quest.negbin.gov.b <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                       parliamentary_office + executive_office + zombie_now + kobo * government, 
                                     data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(quest.negbin.gov.b$model)
quest.negbin.gov.b.se <- coeftest(quest.negbin.gov.b, vcov = vcovCL(quest.negbin.gov.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(quest.negbin.sen.b <- glm.nb(questions ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                       parliamentary_office + executive_office + zombie_now + kobo  * term + government, 
                                     data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(quest.negbin.sen.b$model)
quest.negbin.sen.b.se <- coeftest(quest.negbin.sen.b, vcov = vcovCL(quest.negbin.sen.b, cluster = dat[rows,"candidate_id"]))[,2]

stargazer(list(quest.negbin.a,quest.negbin.gov.a, quest.negbin.b, quest.negbin.sen.b, quest.negbin.gov.b),
          se = list(quest.negbin.a.se,quest.negbin.gov.a.se,quest.negbin.b.se, quest.negbin.sen.b.se, quest.negbin.gov.b.se))


###############################################################################
#                               Table A15                                     #
###############################################################################

summary(ini.negbin.a <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                 parliamentary_office + executive_office + zombie_now + kobo_f  + government, 
                               data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(ini.negbin.a$model)
ini.negbin.a.se <- coeftest(ini.negbin.a, vcov = vcovCL(ini.negbin.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(ini.negbin.gov.a <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + zombie_now + kobo_f  * government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),])) 
rows <- row.names(ini.negbin.gov.a$model)
ini.negbin.gov.a.se <- coeftest(ini.negbin.gov.a, vcov = vcovCL(ini.negbin.gov.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(ini.negbin.b <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                 parliamentary_office + executive_office + zombie_now + kobo  + government, 
                               data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(ini.negbin.b$model)
ini.negbin.b.se <- coeftest(ini.negbin.b, vcov = vcovCL(ini.negbin.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(ini.negbin.sen.b <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + zombie_now + kobo  * term + government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(ini.negbin.sen.b$model)
ini.negbin.sen.b.se <- coeftest(ini.negbin.sen.b, vcov = vcovCL(ini.negbin.sen.b, cluster = dat[rows,"candidate_id"]))[,2]


summary(ini.negbin.gov.b <- glm.nb(num_pmbs_initiated ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                     parliamentary_office + executive_office + zombie_now + kobo * government, 
                                   data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(ini.negbin.gov.b$model)
ini.negbin.gov.b.se <- coeftest(ini.negbin.gov.b, vcov = vcovCL(ini.negbin.gov.b, cluster = dat[rows,"candidate_id"]))[,2]


stargazer(list(ini.negbin.a, ini.negbin.gov.a, ini.negbin.b, ini.negbin.sen.b, ini.negbin.gov.b),
          se = list(ini.negbin.a.se, ini.negbin.gov.a.se, ini.negbin.b.se, ini.negbin.sen.b.se, ini.negbin.gov.b.se))

###############################################################################
#                               Table A16                                     #
###############################################################################

# Operationalization A
summary(co.negbin.a <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                parliamentary_office + executive_office + zombie_now + kobo_f  + government, 
                              data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(co.negbin.a$model)
co.negbin.a.se <- coeftest(co.negbin.a, vcov = vcovCL(co.negbin.a, cluster = dat[rows,"candidate_id"]))[,2]

summary(co.negbin.gov.a <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                    parliamentary_office + executive_office + zombie_now + kobo_f  * government, 
                                  data = dat[which(dat$party_id %in% kobo.parties),])) 
rows <- row.names(co.negbin.gov.a$model)
co.negbin.gov.a.se <- coeftest(co.negbin.gov.a, vcov = vcovCL(co.negbin.gov.a, cluster = dat[rows,"candidate_id"]))[,2]


summary(co.negbin.b <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                parliamentary_office + executive_office + zombie_now + kobo  + government, 
                              data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(co.negbin.b$model)
co.negbin.b.se <- coeftest(co.negbin.b, vcov = vcovCL(co.negbin.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(co.negbin.gov.b <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                    parliamentary_office + executive_office + zombie_now + kobo * government, 
                                  data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(co.negbin.gov.b$model)
co.negbin.gov.b.se <- coeftest(co.negbin.gov.b, vcov = vcovCL(co.negbin.gov.b, cluster = dat[rows,"candidate_id"]))[,2]

summary(co.negbin.sen.b <- glm.nb(num_pmb_cosponsored ~ term + as.factor(party_id) + e2005 + e2009 + e2012+ party_office +
                                    parliamentary_office + executive_office + zombie_now + kobo  * term + government, 
                                  data = dat[which(dat$party_id %in% kobo.parties),]))
rows <- row.names(co.negbin.sen.b$model)
co.negbin.sen.b.se <- coeftest(co.negbin.sen.b, vcov = vcovCL(co.negbin.sen.b, cluster = dat[rows,"candidate_id"]))[,2]


stargazer(list(co.negbin.a, co.negbin.gov.a, co.negbin.b, co.negbin.sen.b, co.negbin.gov.b),
          se = list(co.negbin.a.se, co.negbin.gov.a.se, co.negbin.b.se, co.negbin.sen.b.se, co.negbin.gov.b.se))



###############################################################################
#                                Figure 2a                                    #
###############################################################################

ame.nb <- function(input = dat, dv = c("questions",
                                       "num_pmb_cosponsored","num_pmbs_initiated"),
                   iv = c("kobo","kobo_district"), seed = 1104, nsim = 1000){
  require(sandwich);require(lmtest);require(MASS)
  if(iv == "kobo_district"){dat.tmp <- input[which(input$party_id %in% kobo.parties & input$kobo ==0),]}
  if(iv == "kobo"){dat.tmp <- input[which(input$party_id %in% kobo.parties),]}
  # paste formulas
  form1 <- formula(paste(dv, "~ term + JRP + LDP + YP + e2005 + e2009 + e2012 + party_office + parliamentary_office + executive_office + margin + government +", iv))
  form2 <- formula(paste(dv, "~ JRP + LDP + YP + e2005 + e2009 + e2012 + party_office + parliamentary_office + executive_office + margin + government + term *", iv))
  form3 <- formula(paste(dv, "~ term + JRP + LDP + YP + e2005 + e2009 + e2012 + party_office + parliamentary_office + executive_office + margin + government *", iv))
  # run models
  mod1 <- glm.nb(form1, data = dat.tmp)
  mod2 <- glm.nb(form2, data = dat.tmp)
  mod3 <- glm.nb(form3, data = dat.tmp)
  
  rows1 <- row.names(mod1$model)
  rows2 <- row.names(mod2$model)
  rows3 <- row.names(mod3$model)
  # get variance covariance with clustered variance structure
  mod1.vcov <- vcovCL(mod1, cluster = dat[rows1,"candidate_id"])
  mod2.vcov <- vcovCL(mod2, cluster = dat[rows2,"candidate_id"])
  mod3.vcov <- vcovCL(mod3, cluster = dat[rows3,"candidate_id"])
  # get clustered standard errors
  mod1.se <- sqrt(diag(mod1.vcov))
  mod2.se <- sqrt(diag(mod2.vcov))
  mod3.se <- sqrt(diag(mod3.vcov))
  
  # simulations mod1
  set.seed(seed)
  b.sim1 <- mvrnorm(nsim, coef(mod1), mod1.vcov)
  tmp.df.mod1.0 <- mod1[["model"]][-1]  
  tmp.df.mod1.1 <- mod1[["model"]][-1]
  tmp.df.mod1.0[[iv]] <- 0
  tmp.df.mod1.1[[iv]] <- 1
  
  mod1.pred.0 <- as.matrix(cbind(1,tmp.df.mod1.0)) %*% t(b.sim1)
  mod1.pred.1 <- as.matrix(cbind(1,tmp.df.mod1.1)) %*% t(b.sim1) 
  mod1.pred.0 <- exp(mod1.pred.0)
  mod1.pred.1 <- exp(mod1.pred.1)
  mod1.pred.0.mean <- apply(mod1.pred.0, 2, mean, na.rm = TRUE)
  mod1.pred.1.mean <- apply(mod1.pred.1, 2, mean, na.rm = TRUE)
  mod1.ame.95 <- quantile(mod1.pred.1.mean - mod1.pred.0.mean, probs = c(0.5, 0.025, 0.975), na.rm = TRUE)
  mod1.ame.90 <- quantile(mod1.pred.1.mean - mod1.pred.0.mean, probs = c(0.5, 0.05, 0.95), na.rm = TRUE)
  
  # simulations mod2
  set.seed(seed)
  b.sim2 <- mvrnorm(nsim, coef(mod2), mod2.vcov)
  tmp.df.mod2.0.term.1 <- mod2[["model"]][-1]
  tmp.df.mod2.0.term.2 <- mod2[["model"]][-1]
  tmp.df.mod2.0.term.3 <- mod2[["model"]][-1]
  tmp.df.mod2.1.term.1 <- mod2[["model"]][-1]
  tmp.df.mod2.1.term.2 <- mod2[["model"]][-1]
  tmp.df.mod2.1.term.3 <- mod2[["model"]][-1]
  tmp.df.mod2.0.term.1[[iv]] <- 0 
  tmp.df.mod2.0.term.2[[iv]] <- 0
  tmp.df.mod2.0.term.3[[iv]] <- 0
  tmp.df.mod2.1.term.1[[iv]] <- 1
  tmp.df.mod2.1.term.2[[iv]] <- 1
  tmp.df.mod2.1.term.3[[iv]] <- 1
  
  tmp.df.mod2.0.term.1[["term"]] <- 1 
  tmp.df.mod2.0.term.2[["term"]] <- 2
  tmp.df.mod2.0.term.3[["term"]] <- 3
  tmp.df.mod2.1.term.1[["term"]] <- 1
  tmp.df.mod2.1.term.2[["term"]] <- 2
  tmp.df.mod2.1.term.3[["term"]] <- 3
  
  tmp.df.mod2.0.term.1[[paste0("term:",iv)]] <- 0 
  tmp.df.mod2.0.term.2[[paste0("term:",iv)]] <- 0
  tmp.df.mod2.0.term.3[[paste0("term:",iv)]] <- 0
  tmp.df.mod2.1.term.1[[paste0("term:",iv)]] <- 1
  tmp.df.mod2.1.term.2[[paste0("term:",iv)]] <- 2
  tmp.df.mod2.1.term.3[[paste0("term:",iv)]] <- 3
  
  
  mod2.pred.0.term.1 <- as.matrix(cbind(1,tmp.df.mod2.0.term.1)) %*% t(b.sim2)
  mod2.pred.0.term.2 <- as.matrix(cbind(1,tmp.df.mod2.0.term.2)) %*% t(b.sim2)
  mod2.pred.0.term.3 <- as.matrix(cbind(1,tmp.df.mod2.0.term.3)) %*% t(b.sim2)
  mod2.pred.1.term.1 <- as.matrix(cbind(1,tmp.df.mod2.1.term.1)) %*% t(b.sim2)
  mod2.pred.1.term.2 <- as.matrix(cbind(1,tmp.df.mod2.1.term.2)) %*% t(b.sim2)
  mod2.pred.1.term.3 <- as.matrix(cbind(1,tmp.df.mod2.1.term.3)) %*% t(b.sim2)
  
  mod2.pred.0.term.1 <- exp(mod2.pred.0.term.1) 
  mod2.pred.0.term.2 <- exp(mod2.pred.0.term.2) 
  mod2.pred.0.term.3 <- exp(mod2.pred.0.term.3) 
  mod2.pred.1.term.1 <- exp(mod2.pred.1.term.1)
  mod2.pred.1.term.2 <- exp(mod2.pred.1.term.2)
  mod2.pred.1.term.3 <- exp(mod2.pred.1.term.3)
  
  mod2.pred.0.term.1.mean <- apply(mod2.pred.0.term.1, 2, mean, na.rm = TRUE)
  mod2.pred.0.term.2.mean <- apply(mod2.pred.0.term.2, 2, mean, na.rm = TRUE)
  mod2.pred.0.term.3.mean <- apply(mod2.pred.0.term.3, 2, mean, na.rm = TRUE)
  mod2.pred.1.term.1.mean <- apply(mod2.pred.1.term.1, 2, mean, na.rm = TRUE)
  mod2.pred.1.term.2.mean <- apply(mod2.pred.1.term.2, 2, mean, na.rm = TRUE)
  mod2.pred.1.term.3.mean <- apply(mod2.pred.1.term.3, 2, mean, na.rm = TRUE)
  
  mod2.term.1.ame.95 <- quantile(mod2.pred.1.term.1.mean - mod2.pred.0.term.1.mean, probs = c(0.5, 0.025, 0.975), na.rm = TRUE)
  mod2.term.2.ame.95 <- quantile(mod2.pred.1.term.2.mean - mod2.pred.0.term.2.mean, probs = c(0.5, 0.025, 0.975), na.rm = TRUE)
  mod2.term.3.ame.95 <- quantile(mod2.pred.1.term.3.mean - mod2.pred.0.term.3.mean, probs = c(0.5, 0.025, 0.975), na.rm = TRUE)
  
  mod2.term.1.ame.90 <- quantile(mod2.pred.1.term.1.mean - mod2.pred.0.term.1.mean, probs = c(0.5, 0.05, 0.95), na.rm = TRUE)
  mod2.term.2.ame.90 <- quantile(mod2.pred.1.term.2.mean - mod2.pred.0.term.2.mean, probs = c(0.5, 0.05, 0.95), na.rm = TRUE)
  mod2.term.3.ame.90 <- quantile(mod2.pred.1.term.3.mean - mod2.pred.0.term.3.mean, probs = c(0.5, 0.05, 0.95), na.rm = TRUE)
  
  # simulations mod3
  set.seed(seed)
  b.sim3 <- mvrnorm(nsim, coef(mod3), mod3.vcov)
  tmp.df.mod3.0.gov.0 <- mod3[["model"]][-1]
  tmp.df.mod3.0.gov.1 <- mod3[["model"]][-1]
  tmp.df.mod3.1.gov.0 <- mod3[["model"]][-1]
  tmp.df.mod3.1.gov.1 <- mod3[["model"]][-1]
  
  tmp.df.mod3.0.gov.0[[iv]] <- 0
  tmp.df.mod3.0.gov.1[[iv]] <- 0
  tmp.df.mod3.1.gov.0[[iv]] <- 1
  tmp.df.mod3.1.gov.1[[iv]] <- 1
  
  tmp.df.mod3.0.gov.0[["government"]] <- 0
  tmp.df.mod3.0.gov.1[["government"]] <- 1
  tmp.df.mod3.1.gov.0[["government"]] <- 0
  tmp.df.mod3.1.gov.1[["government"]] <- 1
  
  tmp.df.mod3.0.gov.0[[paste0("government:",iv)]] <- 0
  tmp.df.mod3.0.gov.1[[paste0("government:",iv)]] <- 0
  tmp.df.mod3.1.gov.0[[paste0("government:",iv)]] <- 0
  tmp.df.mod3.1.gov.1[[paste0("government:",iv)]] <- 1
  
  mod3.pred.0.gov.0 <- as.matrix(cbind(1,tmp.df.mod3.0.gov.0)) %*% t(b.sim3)
  mod3.pred.0.gov.1 <- as.matrix(cbind(1,tmp.df.mod3.0.gov.1)) %*% t(b.sim3)
  mod3.pred.1.gov.0 <- as.matrix(cbind(1,tmp.df.mod3.1.gov.0)) %*% t(b.sim3)
  mod3.pred.1.gov.1 <- as.matrix(cbind(1,tmp.df.mod3.1.gov.1)) %*% t(b.sim3)
  
  mod3.pred.0.gov.0 <- exp(mod3.pred.0.gov.0)
  mod3.pred.0.gov.1 <- exp(mod3.pred.0.gov.1)
  mod3.pred.1.gov.0 <- exp(mod3.pred.1.gov.0)
  mod3.pred.1.gov.1 <- exp(mod3.pred.1.gov.1)
  
  mod3.pred.0.gov.0.mean <- apply(mod3.pred.0.gov.0, 2, mean, na.rm = TRUE)
  mod3.pred.0.gov.1.mean <- apply(mod3.pred.0.gov.1, 2, mean, na.rm = TRUE)
  mod3.pred.1.gov.0.mean <- apply(mod3.pred.1.gov.0, 2, mean, na.rm = TRUE)
  mod3.pred.1.gov.1.mean <- apply(mod3.pred.1.gov.1, 2, mean, na.rm = TRUE)
  
  mod3.gov.0.ame.95 <- quantile(mod3.pred.1.gov.0.mean - mod3.pred.0.gov.0.mean, probs = c(0.5, 0.025, 0.975), na.rm = TRUE)
  mod3.gov.1.ame.95 <- quantile(mod3.pred.1.gov.1.mean - mod3.pred.0.gov.1.mean, probs = c(0.5, 0.025, 0.975), na.rm = TRUE)
  
  mod3.gov.0.ame.90 <- quantile(mod3.pred.1.gov.0.mean - mod3.pred.0.gov.0.mean, probs = c(0.5, 0.05, 0.95), na.rm = TRUE)
  mod3.gov.1.ame.90 <- quantile(mod3.pred.1.gov.1.mean - mod3.pred.0.gov.1.mean, probs = c(0.5, 0.05, 0.95), na.rm = TRUE)
  
  return(list(main_model = mod1, senior_model = mod2, gov_model = mod3, 
              main_vcov = mod1.vcov, senior_vcov = mod2.vcov, gov_model = mod3.vcov, 
              main_se = mod1.se, senior_se = mod2.se, gov_se = mod3.se,
              ame_main90 = mod1.ame.90, ame_main95 = mod1.ame.95,
              ame_senior1_90 = mod2.term.1.ame.90, ame_senior2_90 = mod2.term.2.ame.90,
              ame_senior3_90 = mod2.term.3.ame.90, ame_senior1_95 = mod2.term.1.ame.95, 
              ame_senior2_95 = mod2.term.2.ame.95, ame_senior3_95 = mod2.term.3.ame.95,
              ame_gov0_90 = mod3.gov.0.ame.90, ame_gov1_90 = mod3.gov.1.ame.90,
              ame_gov0_95 = mod3.gov.0.ame.95, ame_gov1_95 = mod3.gov.1.ame.95))
}


ame.quest.robust <- ame.nb(input = dat[dat$kobo == 0,], dv = "questions", iv = "kobo_district", 1104, 1000)
ame.ini.robust <- ame.nb(input = dat[dat$kobo == 0,], dv = "num_pmbs_initiated", iv = "kobo_district", 1104, 1000)
ame.co.robust <- ame.nb(input = dat[dat$kobo == 0,], dv = "num_pmb_cosponsored", iv = "kobo_district", 1104, 1000)

ame.quest <- ame.nb(input = dat, dv = "questions", iv = "kobo", 1104, 1000)
ame.ini <- ame.nb(input = dat, dv = "num_pmbs_initiated", iv = "kobo", 1104, 1000)
ame.co <- ame.nb(input = dat, dv = "num_pmb_cosponsored", iv = "kobo", 1104, 1000)

ggData <- as.data.frame(rbind(cbind(iv = "kobo_district", dv = "Written   \nQuestions ", ci = 90, ame = ame.quest.robust$ame_main90[1], upper = ame.quest.robust$ame_main90[3],lower = ame.quest.robust$ame_main90[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nInitiated   ", ci = 90, ame = ame.ini.robust$ame_main90[1], upper = ame.ini.robust$ame_main90[3],lower = ame.ini.robust$ame_main90[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nCosponsored", ci = 90, ame = ame.co.robust$ame_main90[1], upper = ame.co.robust$ame_main90[3],lower = ame.co.robust$ame_main90[2]),
  cbind(iv = "kobo_district", dv = "Written   \nQuestions ", ci = 95, ame = ame.quest.robust$ame_main95[1], upper = ame.quest.robust$ame_main95[3],lower = ame.quest.robust$ame_main95[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nInitiated   ", ci = 95, ame = ame.ini.robust$ame_main95[1], upper = ame.ini.robust$ame_main95[3],lower = ame.ini.robust$ame_main95[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nCosponsored", ci = 95, ame = ame.co.robust$ame_main95[1], upper = ame.co.robust$ame_main95[3],lower = ame.co.robust$ame_main95[2]),
  cbind(iv = "kobo", dv = "Written   \nQuestions ", ci = 90, ame = ame.quest$ame_main90[1], upper = ame.quest$ame_main90[3],lower = ame.quest$ame_main90[2]),
  cbind(iv = "kobo", dv = "PMBs     \nInitiated   ", ci = 90, ame = ame.ini$ame_main90[1], upper = ame.ini$ame_main90[3],lower = ame.ini$ame_main90[2]),
  cbind(iv = "kobo", dv = "PMBs     \nCosponsored", ci = 90, ame = ame.co$ame_main90[1], upper = ame.co$ame_main90[3],lower = ame.co$ame_main90[2]),
  cbind(iv = "kobo", dv = "Written   \nQuestions ", ci = 95, ame = ame.quest$ame_main95[1], upper = ame.quest$ame_main95[3],lower = ame.quest$ame_main95[2]),
  cbind(iv = "kobo", dv = "PMBs     \nInitiated   ", ci = 95, ame = ame.ini$ame_main95[1], upper = ame.ini$ame_main95[3],lower = ame.ini$ame_main95[2]),
  cbind(iv = "kobo", dv = "PMBs     \nCosponsored", ci = 95, ame = ame.co$ame_main95[1], upper = ame.co$ame_main95[3],lower = ame.co$ame_main95[2])))

ggData[,4] <- as.numeric(as.character(ggData[,4]))
ggData[,5] <- as.numeric(as.character(ggData[,5]))
ggData[,6] <- as.numeric(as.character(ggData[,6]))
ggData$sig <- as.numeric(0)
ggData$sig[(ggData$upper> 0 & ggData$lower > 0 )|(ggData$upper< 0 & ggData$lower < 0 )  ] <- 1


yticks <- c(1.8,1.5,1.2)
xticks <- c(-3,-2,-1,0)
pdf("C:/Users/rehmerjo/Dropbox/HERTIE/paper_kobo/R&R/tex/kobo_main.pdf")
par(oma = c(3,1,1,1),
    mai = c(3,1,1,1),
    mar =  c(3, 5.5, 1, 1) + 0.1,
    mgp = c(3, 1, 1), 
    xpd = NA,
    bty = 'n')
plot(ggData$ame[ggData$ci == 95 & ggData$iv == "kobo"], yticks, 
     ylab = "", xlab = "", pch = 19, xlim = c(-3,0), ylim = c(1,2), 
     cex.lab = 1.6, frame = FALSE, xaxt='n', yaxt = 'n')
segments( ggData$lower[ggData$iv == "kobo" & ggData$ci == 95], yticks,
         ggData$upper[ggData$iv == "kobo" & ggData$ci == 95] , yticks, lty = 2   )
segments( ggData$lower[ggData$iv == "kobo" & ggData$ci == 90], yticks,
          ggData$upper[ggData$iv == "kobo" & ggData$ci == 90] , yticks, lty = 1   )
text(x = rep(-3.3,3), y = yticks, labels = c("PQs","PMBs init.","PMBs cosp."))
lines(x = rep(0,3), y = c(1.1,1.5,1.9), col = "red", lty = "dashed")
text(par("usr")[3], x = xticks, labels = xticks, 
     pos = 1, xpd = TRUE, cex = 1.25)
axis(side = 1, at = -1.5,labels = "Average Marginal Effects", col = "white", cex.axis = 1.5, line = 1)
dev.off()

###############################################################################
#                                Figure 2b                                    #
###############################################################################

yticks <- c(1.8,1.5,1.2)
xticks <- c(-2,0,2,4)
#pdf("C:/Users/rehmerjo/Dropbox/HERTIE/paper_kobo/R&R/tex/kobo_district_main.pdf")
par(oma = c(3,1,1,1),
    mai = c(3,1,1,1),
    mar =  c(3, 5.5, 1, 1) + 0.1,
    mgp = c(3, 1, 1), 
    xpd = NA,
    bty = 'n')
plot(ggData$ame[ggData$ci == 95 & ggData$iv == "kobo_district"], yticks, 
     ylab = "", xlab = "", pch = 19, xlim = c(-2,4), ylim = c(1,2), 
     cex.lab = 1.6, frame = FALSE, xaxt='n', yaxt = 'n')
segments( ggData$lower[ggData$iv == "kobo_district" & ggData$ci == 95], yticks,
          ggData$upper[ggData$iv == "kobo_district" & ggData$ci == 95] , yticks, lty = 2   )
segments( ggData$lower[ggData$iv == "kobo_district" & ggData$ci == 90], yticks,
          ggData$upper[ggData$iv == "kobo_district" & ggData$ci == 90] , yticks, lty = 1   )
lines(x = rep(0,3), y = c(1.1,1.5,1.9), col = "red", lty = "dashed")
text(par("usr")[3], x = xticks, labels = xticks, 
     pos = 1, xpd = TRUE, cex = 1.25)
axis(side = 1, at = 1,labels = "Average Marginal Effects", col = "white", cex.axis = 1.5, line = 1)
#dev.off()


###############################################################################
#                                Figure 3                                     #
###############################################################################

ggData <- as.data.frame(rbind(cbind(iv = "kobo_district", dv = "Written   \nQuestions ", term = 1, ci = 90, ame = ame.quest.robust$ame_senior1_90[1], upper = ame.quest.robust$ame_senior1_90[3],lower = ame.quest.robust$ame_senior1_90[2]),
  cbind(iv = "kobo_district", dv = "Written   \nQuestions ", term = 2, ci = 90, ame = ame.quest.robust$ame_senior2_90[1], upper = ame.quest.robust$ame_senior2_90[3],lower = ame.quest.robust$ame_senior2_90[2]),
  cbind(iv = "kobo_district", dv = "Written   \nQuestions ", term = 3, ci = 90, ame = ame.quest.robust$ame_senior3_90[1], upper = ame.quest.robust$ame_senior3_90[3],lower = ame.quest.robust$ame_senior3_90[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nInitiated   ", term = 1, ci = 90, ame = ame.ini.robust$ame_senior1_90[1], upper = ame.ini.robust$ame_senior1_90[3],lower = ame.ini.robust$ame_senior1_90[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nInitiated   ", term = 2, ci = 90, ame = ame.ini.robust$ame_senior2_90[1], upper = ame.ini.robust$ame_senior2_90[3],lower = ame.ini.robust$ame_senior2_90[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nInitiated   ", term = 3, ci = 90, ame = ame.ini.robust$ame_senior3_90[1], upper = ame.ini.robust$ame_senior3_90[3],lower = ame.ini.robust$ame_senior3_90[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nCosponsored", term = 1, ci = 90, ame = ame.co.robust$ame_senior1_90[1], upper = ame.co.robust$ame_senior1_90[3],lower = ame.co.robust$ame_senior1_90[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nCosponsored", term = 2, ci = 90, ame = ame.co.robust$ame_senior2_90[1], upper = ame.co.robust$ame_senior2_90[3],lower = ame.co.robust$ame_senior2_90[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nCosponsored", term = 3, ci = 90, ame = ame.co.robust$ame_senior3_90[1], upper = ame.co.robust$ame_senior3_90[3],lower = ame.co.robust$ame_senior3_90[2]),
  cbind(iv = "kobo_district", dv = "Written   \nQuestions ", term = 1, ci = 95, ame = ame.quest.robust$ame_senior1_95[1], upper = ame.quest.robust$ame_senior1_95[3],lower = ame.quest.robust$ame_senior1_95[2]),
  cbind(iv = "kobo_district", dv = "Written   \nQuestions ", term = 2, ci = 95, ame = ame.quest.robust$ame_senior2_95[1], upper = ame.quest.robust$ame_senior2_95[3],lower = ame.quest.robust$ame_senior2_95[2]),
  cbind(iv = "kobo_district", dv = "Written   \nQuestions ", term = 3, ci = 95, ame = ame.quest.robust$ame_senior3_95[1], upper = ame.quest.robust$ame_senior3_95[3],lower = ame.quest.robust$ame_senior3_95[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nInitiated   ", term = 1, ci = 95, ame = ame.ini.robust$ame_senior1_95[1], upper = ame.ini.robust$ame_senior1_95[3],lower = ame.ini.robust$ame_senior1_95[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nInitiated   ", term = 2, ci = 95, ame = ame.ini.robust$ame_senior2_95[1], upper = ame.ini.robust$ame_senior2_95[3],lower = ame.ini.robust$ame_senior2_95[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nInitiated   ", term = 3, ci = 95, ame = ame.ini.robust$ame_senior3_95[1], upper = ame.ini.robust$ame_senior3_95[3],lower = ame.ini.robust$ame_senior3_95[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nCosponsored", term = 1, ci = 95, ame = ame.co.robust$ame_senior1_95[1], upper = ame.co.robust$ame_senior1_95[3],lower = ame.co.robust$ame_senior1_95[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nCosponsored", term = 2, ci = 95, ame = ame.co.robust$ame_senior2_95[1], upper = ame.co.robust$ame_senior2_95[3],lower = ame.co.robust$ame_senior2_95[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nCosponsored", term = 3, ci = 95, ame = ame.co.robust$ame_senior3_95[1], upper = ame.co.robust$ame_senior3_95[3],lower = ame.co.robust$ame_senior3_95[2]),
  cbind(iv = "kobo", dv = "Written   \nQuestions ", term = 1, ci = 90, ame = ame.quest$ame_senior1_90[1], upper = ame.quest$ame_senior1_90[3],lower = ame.quest$ame_senior1_90[2]),
  cbind(iv = "kobo", dv = "Written   \nQuestions ", term = 2, ci = 90, ame = ame.quest$ame_senior2_90[1], upper = ame.quest$ame_senior2_90[3],lower = ame.quest$ame_senior2_90[2]),
  cbind(iv = "kobo", dv = "Written   \nQuestions ", term = 3, ci = 90, ame = ame.quest$ame_senior3_90[1], upper = ame.quest$ame_senior3_90[3],lower = ame.quest$ame_senior3_90[2]),
  cbind(iv = "kobo", dv = "PMBs     \nInitiated   ", term = 1, ci = 90, ame = ame.ini$ame_senior1_90[1], upper = ame.ini$ame_senior1_90[3],lower = ame.ini$ame_senior1_90[2]),
  cbind(iv = "kobo", dv = "PMBs     \nInitiated   ", term = 2, ci = 90, ame = ame.ini$ame_senior2_90[1], upper = ame.ini$ame_senior2_90[3],lower = ame.ini$ame_senior2_90[2]),
  cbind(iv = "kobo", dv = "PMBs     \nInitiated   ", term = 3, ci = 90, ame = ame.ini$ame_senior3_90[1], upper = ame.ini$ame_senior3_90[3],lower = ame.ini$ame_senior3_90[2]),
  cbind(iv = "kobo", dv = "PMBs     \nCosponsored", term = 1, ci = 90, ame = ame.co$ame_senior1_90[1], upper = ame.co$ame_senior1_90[3],lower = ame.co$ame_senior1_90[2]),
  cbind(iv = "kobo", dv = "PMBs     \nCosponsored", term = 2, ci = 90, ame = ame.co$ame_senior2_90[1], upper = ame.co$ame_senior2_90[3],lower = ame.co$ame_senior2_90[2]),
  cbind(iv = "kobo", dv = "PMBs     \nCosponsored", term = 3, ci = 90, ame = ame.co$ame_senior3_90[1], upper = ame.co$ame_senior3_90[3],lower = ame.co$ame_senior3_90[2]),
  cbind(iv = "kobo", dv = "Written   \nQuestions ", term = 1, ci = 95, ame = ame.quest$ame_senior1_95[1], upper = ame.quest$ame_senior1_95[3],lower = ame.quest$ame_senior1_95[2]),
  cbind(iv = "kobo", dv = "Written   \nQuestions ", term = 2, ci = 95, ame = ame.quest$ame_senior2_95[1], upper = ame.quest$ame_senior2_95[3],lower = ame.quest$ame_senior2_95[2]),
  cbind(iv = "kobo", dv = "Written   \nQuestions ", term = 3, ci = 95, ame = ame.quest$ame_senior3_95[1], upper = ame.quest$ame_senior3_95[3],lower = ame.quest$ame_senior3_95[2]),
  cbind(iv = "kobo", dv = "PMBs     \nInitiated   ", term = 1, ci = 95, ame = ame.ini$ame_senior1_95[1], upper = ame.ini$ame_senior1_95[3],lower = ame.ini$ame_senior1_95[2]),
  cbind(iv = "kobo", dv = "PMBs     \nInitiated   ", term = 2, ci = 95, ame = ame.ini$ame_senior2_95[1], upper = ame.ini$ame_senior2_95[3],lower = ame.ini$ame_senior2_95[2]),
  cbind(iv = "kobo", dv = "PMBs     \nInitiated   ", term = 3, ci = 95, ame = ame.ini$ame_senior3_95[1], upper = ame.ini$ame_senior3_95[3],lower = ame.ini$ame_senior3_95[2]),
  cbind(iv = "kobo", dv = "PMBs     \nCosponsored", term = 1, ci = 95, ame = ame.co$ame_senior1_95[1], upper = ame.co$ame_senior1_95[3],lower = ame.co$ame_senior1_95[2]),
  cbind(iv = "kobo", dv = "PMBs     \nCosponsored", term = 2, ci = 95, ame = ame.co$ame_senior2_95[1], upper = ame.co$ame_senior2_95[3],lower = ame.co$ame_senior2_95[2]),
  cbind(iv = "kobo", dv = "PMBs     \nCosponsored", term = 3, ci = 95, ame = ame.co$ame_senior3_95[1], upper = ame.co$ame_senior3_95[3],lower = ame.co$ame_senior3_95[2])))

ggData[,7] <- as.numeric(as.character(ggData[,7]))
ggData[,5] <- as.numeric(as.character(ggData[,5]))
ggData[,6] <- as.numeric(as.character(ggData[,6]))
ggData$sig <- as.numeric(0)
ggData$sig[(ggData$upper> 0 & ggData$lower > 0 )|(ggData$upper< 0 & ggData$lower < 0 )  ] <- 1
ggData$term_f <- ""
ggData$term_f[ggData$term == "1"] <- "1st Term"
ggData$term_f[ggData$term == "2"] <- "2nd Term"
ggData$term_f[ggData$term == "3"] <- "3rd Term"

###############################################################################
#                                Figure 3a                                     #
###############################################################################


yticks <- c(1.8,1.5,1.2)
xticks <- c(-5,-3,-1,0)
pdf("C:/Users/rehmerjo/Dropbox/HERTIE/paper_kobo/R&R/tex/kobo_senior.pdf",width = 12, height = 6)
par(mfrow = c(1,3),
    oma = c(3,1,1,1),
    mai = c(1, 0.1, 0.1, 0.1),
    mar =  c(2, 5, 0, 1) + 0.1,
    mgp = c(3, 1, 1), 
    xpd = NA,
    bty = 'n')
plot(ggData$ame[ggData$ci == 95 & ggData$iv == "kobo" & ggData$term == 1], yticks, 
     ylab = "", xlab = "", pch = 19, xlim = c(-5,1), ylim = c(1.1,1.9), 
     cex.lab = 1.6, frame = FALSE, xaxt='n', yaxt = 'n', main = "1st Term", cex.main = 1.75)
segments( ggData$lower[ggData$iv == "kobo" & ggData$ci == 95 & ggData$term == 1], yticks,
          ggData$upper[ggData$iv == "kobo" & ggData$ci == 95 & ggData$term == 1] , yticks, lty = 2   )
segments( ggData$lower[ggData$iv == "kobo" & ggData$ci == 90 & ggData$term == 1], yticks,
          ggData$upper[ggData$iv == "kobo" & ggData$ci == 90 & ggData$term == 1] , yticks, lty = 1   )
text(x = rep(-5.4,3), y = yticks, labels = c("PQs","PMBs init.","PMBs cosp."), cex = 1.5)
lines(x = rep(0,3), y = c(1.1,1.5,1.9), col = "red", lty = "dashed")
text(par("usr")[3], x = xticks, labels = xticks, 
     pos = 3, xpd = TRUE, cex = 1.25)
axis(side = 1, at = -2,labels = "Average Marginal Effects", col = "white", cex.axis = 1.75, line = 0)


plot(ggData$ame[ggData$ci == 95 & ggData$iv == "kobo" & ggData$term == 2], yticks, 
     ylab = "", xlab = "", pch = 19, xlim = c(-5,1), ylim = c(1.1,1.9), 
     cex.lab = 1.6, frame = FALSE, xaxt='n', yaxt = 'n' , main = "2nd Term", cex.main = 1.75)
segments( ggData$lower[ggData$iv == "kobo" & ggData$ci == 95 & ggData$term == 2], yticks,
          ggData$upper[ggData$iv == "kobo" & ggData$ci == 95 & ggData$term == 2] , yticks, lty = 2   )
segments( ggData$lower[ggData$iv == "kobo" & ggData$ci == 90 & ggData$term == 2], yticks,
          ggData$upper[ggData$iv == "kobo" & ggData$ci == 90 & ggData$term == 2] , yticks, lty = 1   )
#text(x = rep(-5.3,3), y = yticks, labels = c("PQs","PMBs init.","PMBs cosp."))
lines(x = rep(0,3), y = c(1.1,1.5,1.9), col = "red", lty = "dashed")
text(par("usr")[3], x = xticks, labels = xticks, 
     pos = 3, xpd = TRUE, cex = 1.25)
axis(side = 1, at = -2,labels = "Average Marginal Effects", col = "white", cex.axis = 1.75, line = 0)


plot(ggData$ame[ggData$ci == 95 & ggData$iv == "kobo" & ggData$term == 3], yticks, 
     ylab = "", xlab = "", pch = 19, xlim = c(-5,1), ylim = c(1.1,1.9), 
     cex.lab = 1.6, frame = FALSE, xaxt='n', yaxt = 'n', main = "3rd Term", cex.main = 1.75)
segments( ggData$lower[ggData$iv == "kobo" & ggData$ci == 95 & ggData$term == 3], yticks,
          ggData$upper[ggData$iv == "kobo" & ggData$ci == 95 & ggData$term == 3] , yticks, lty = 2   )
segments( ggData$lower[ggData$iv == "kobo" & ggData$ci == 90 & ggData$term == 3], yticks,
          ggData$upper[ggData$iv == "kobo" & ggData$ci == 90 & ggData$term == 3] , yticks, lty = 1   )
#text(x = rep(-5.3,3), y = yticks, labels = c("PQs","PMBs init.","PMBs cosp."))
lines(x = rep(0,3), y = c(1.1,1.5,1.9), col = "red", lty = "dashed")
text(par("usr")[3], x = xticks, labels = xticks, 
     pos = 3, xpd = TRUE, cex = 1.25)
axis(side = 1, at = -2,labels = "Average Marginal Effects", col = "white", cex.axis = 1.75, line = 0)
dev.off()

###############################################################################
#                                Figure 3b                                    #
###############################################################################


yticks <- c(1.8,1.5,1.2)
xticks <- c(-2,0,2,4,6,8,10)
pdf("C:/Users/rehmerjo/Dropbox/HERTIE/paper_kobo/R&R/tex/kobo_district_senior.pdf",width = 12, height = 6)
par(mfrow = c(1,3),
    oma = c(3,1,1,1),
    mai = c(1, 0.1, 0.1, 0.1),
    mar =  c(2, 5, 0, 1) + 0.1,
    mgp = c(3, 1, 1), 
    xpd = NA,
    bty = 'n')
plot(ggData$ame[ggData$ci == 95 & ggData$iv == "kobo_district" & ggData$term == 1], yticks, 
     ylab = "", xlab = "", pch = 19, xlim = c(-5,10), ylim = c(1.1,1.9), 
     cex.lab = 1.6, frame = FALSE, xaxt='n', yaxt = 'n', main = "1st Term", cex.main = 1.75)
segments( ggData$lower[ggData$iv == "kobo_district" & ggData$ci == 95 & ggData$term == 1], yticks,
          ggData$upper[ggData$iv == "kobo_district" & ggData$ci == 95 & ggData$term == 1] , yticks, lty = 2   )
segments( ggData$lower[ggData$iv == "kobo_district" & ggData$ci == 90 & ggData$term == 1], yticks,
          ggData$upper[ggData$iv == "kobo_district" & ggData$ci == 90 & ggData$term == 1] , yticks, lty = 1   )
text(x = rep(-5.4,3), y = yticks, labels = c("PQs","PMBs init.","PMBs cosp."), cex = 1.5)
lines(x = rep(0,3), y = c(1.1,1.5,1.9), col = "red", lty = "dashed")
text(par("usr")[3], x = xticks, labels = xticks, 
     pos = 3, xpd = TRUE, cex = 1.25)
axis(side = 1, at = 4,labels = "Average Marginal Effects", col = "white", cex.axis = 1.75, line = 0)


plot(ggData$ame[ggData$ci == 95 & ggData$iv == "kobo_district" & ggData$term == 2], yticks, 
     ylab = "", xlab = "", pch = 19, xlim = c(-5,10), ylim = c(1.1,1.9), 
     cex.lab = 1.6, frame = FALSE, xaxt='n', yaxt = 'n' , main = "2nd Term", cex.main = 1.75)
segments( ggData$lower[ggData$iv == "kobo_district" & ggData$ci == 95 & ggData$term == 2], yticks,
          ggData$upper[ggData$iv == "kobo_district" & ggData$ci == 95 & ggData$term == 2] , yticks, lty = 2   )
segments( ggData$lower[ggData$iv == "kobo_district" & ggData$ci == 90 & ggData$term == 2], yticks,
          ggData$upper[ggData$iv == "kobo_district" & ggData$ci == 90 & ggData$term == 2] , yticks, lty = 1   )
#text(x = rep(-5.3,3), y = yticks, labels = c("PQs","PMBs init.","PMBs cosp."))
lines(x = rep(0,3), y = c(1.1,1.5,1.9), col = "red", lty = "dashed")
text(par("usr")[3], x = xticks, labels = xticks, 
     pos = 3, xpd = TRUE, cex = 1.25)
axis(side = 1, at = 4,labels = "Average Marginal Effects", col = "white", cex.axis = 1.75, line = 0)


plot(ggData$ame[ggData$ci == 95 & ggData$iv == "kobo_district" & ggData$term == 3], yticks, 
     ylab = "", xlab = "", pch = 19, xlim = c(-5,10), ylim = c(1.1,1.9), 
     cex.lab = 1.6, frame = FALSE, xaxt='n', yaxt = 'n', main = "3rd Term", cex.main = 1.75)
segments( ggData$lower[ggData$iv == "kobo_district" & ggData$ci == 95 & ggData$term == 3], yticks,
          ggData$upper[ggData$iv == "kobo_district" & ggData$ci == 95 & ggData$term == 3] , yticks, lty = 2   )
segments( ggData$lower[ggData$iv == "kobo_district" & ggData$ci == 90 & ggData$term == 3], yticks,
          ggData$upper[ggData$iv == "kobo_district" & ggData$ci == 90 & ggData$term == 3] , yticks, lty = 1   )
#text(x = rep(-5.3,3), y = yticks, labels = c("PQs","PMBs init.","PMBs cosp."))
lines(x = rep(0,3), y = c(1.1,1.5,1.9), col = "red", lty = "dashed")
text(par("usr")[3], x = xticks, labels = xticks, 
     pos = 3, xpd = TRUE, cex = 1.25)
axis(side = 1, at = 4,labels = "Average Marginal Effects", col = "white", cex.axis = 1.75, line = 0)
dev.off()


###############################################################################
#                                Figure 4                                     #
###############################################################################

ggData <- as.data.frame(rbind(cbind(iv = "kobo_district", dv = "Written   \nQuestions ", ci = 90, gov = 0,ame = ame.quest.robust$ame_gov0_90[1], upper = ame.quest.robust$ame_gov0_90[3],lower = ame.quest.robust$ame_gov0_90[2]),
  cbind(iv = "kobo_district", dv = "Written   \nQuestions ", ci = 90, gov = 1,ame = ame.quest.robust$ame_gov1_90[1], upper = ame.quest.robust$ame_gov1_90[3],lower = ame.quest.robust$ame_gov1_90[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nInitiated   ", ci = 90, gov = 0,ame = ame.ini.robust$ame_gov0_90[1], upper = ame.ini.robust$ame_gov0_90[3],lower = ame.ini.robust$ame_gov0_90[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nInitiated   ", ci = 90, gov = 1,ame = ame.ini.robust$ame_gov1_90[1], upper = ame.ini.robust$ame_gov1_90[3],lower = ame.ini.robust$ame_gov1_90[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nCosponsored", ci = 90, gov = 0,ame = ame.co.robust$ame_gov0_90[1], upper = ame.co.robust$ame_gov0_90[3],lower = ame.co.robust$ame_gov0_90[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nCosponsored", ci = 90, gov = 1,ame = ame.co.robust$ame_gov1_90[1], upper = ame.co.robust$ame_gov1_90[3],lower = ame.co.robust$ame_gov1_90[2]),
  cbind(iv = "kobo_district", dv = "Written   \nQuestions ", ci = 95, gov = 0,ame = ame.quest.robust$ame_gov0_95[1], upper = ame.quest.robust$ame_gov0_95[3],lower = ame.quest.robust$ame_gov0_95[2]),
  cbind(iv = "kobo_district", dv = "Written   \nQuestions ", ci = 95, gov = 1,ame = ame.quest.robust$ame_gov1_95[1], upper = ame.quest.robust$ame_gov1_95[3],lower = ame.quest.robust$ame_gov1_95[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nInitiated   ", ci = 95, gov = 0,ame = ame.ini.robust$ame_gov0_95[1], upper = ame.ini.robust$ame_gov0_95[3],lower = ame.ini.robust$ame_gov0_95[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nInitiated   ", ci = 95, gov = 1,ame = ame.ini.robust$ame_gov1_95[1], upper = ame.ini.robust$ame_gov1_95[3],lower = ame.ini.robust$ame_gov1_95[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nCosponsored", ci = 95, gov = 0,ame = ame.co.robust$ame_gov0_95[1], upper = ame.co.robust$ame_gov0_95[3],lower = ame.co.robust$ame_gov0_95[2]),
  cbind(iv = "kobo_district", dv = "PMBs     \nCosponsored", ci = 95, gov = 1,ame = ame.co.robust$ame_gov1_95[1], upper = ame.co.robust$ame_gov1_95[3],lower = ame.co.robust$ame_gov1_95[2]),
  cbind(iv = "kobo", dv = "Written   \nQuestions ", ci = 90, gov = 1,ame = ame.quest$ame_gov1_90[1], upper = ame.quest$ame_gov1_90[3],lower = ame.quest$ame_gov1_90[2]),
  cbind(iv = "kobo", dv = "Written   \nQuestions ", ci = 90, gov = 0,ame = ame.quest$ame_gov0_90[1], upper = ame.quest$ame_gov0_90[3],lower = ame.quest$ame_gov0_90[2]),
  cbind(iv = "kobo", dv = "PMBs     \nInitiated   ", ci = 90, gov = 1,ame = ame.ini$ame_gov1_90[1], upper = ame.ini$ame_gov1_90[3],lower = ame.ini$ame_gov1_90[2]),
  cbind(iv = "kobo", dv = "PMBs     \nInitiated   ", ci = 90, gov = 0,ame = ame.ini$ame_gov0_90[1], upper = ame.ini$ame_gov0_90[3],lower = ame.ini$ame_gov0_90[2]),
  cbind(iv = "kobo", dv = "PMBs     \nCosponsored", ci = 90, gov = 1,ame = ame.co$ame_gov1_90[1], upper = ame.co$ame_gov1_90[3],lower = ame.co$ame_gov1_90[2]),
  cbind(iv = "kobo", dv = "PMBs     \nCosponsored", ci = 90, gov = 0,ame = ame.co$ame_gov0_90[1], upper = ame.co$ame_gov0_90[3],lower = ame.co$ame_gov0_90[2]),
  cbind(iv = "kobo", dv = "Written   \nQuestions ", ci = 95, gov = 1,ame = ame.quest$ame_gov1_95[1], upper = ame.quest$ame_gov1_95[3],lower = ame.quest$ame_gov1_95[2]),
  cbind(iv = "kobo", dv = "Written   \nQuestions ", ci = 95, gov = 0,ame = ame.quest$ame_gov0_95[1], upper = ame.quest$ame_gov0_95[3],lower = ame.quest$ame_gov0_95[2]),
  cbind(iv = "kobo", dv = "PMBs     \nInitiated   ", ci = 95, gov = 1,ame = ame.ini$ame_gov1_95[1], upper = ame.ini$ame_gov1_95[3],lower = ame.ini$ame_gov1_95[2]),
  cbind(iv = "kobo", dv = "PMBs     \nInitiated   ", ci = 95, gov = 0, ame = ame.ini$ame_gov0_95[1], upper = ame.ini$ame_gov0_95[3],lower = ame.ini$ame_gov0_95[2]),
  cbind(iv = "kobo", dv = "PMBs     \nCosponsored", ci = 95, gov = 1,ame = ame.co$ame_gov1_95[1], upper = ame.co$ame_gov1_95[3],lower = ame.co$ame_gov1_95[2]),
  cbind(iv = "kobo", dv = "PMBs     \nCosponsored", ci = 95, gov = 0, ame = ame.co$ame_gov0_95[1], upper = ame.co$ame_gov0_95[3],lower = ame.co$ame_gov0_95[2]))) 

ggData[,7] <- as.numeric(as.character(ggData[,7]))
ggData[,5] <- as.numeric(as.character(ggData[,5]))
ggData[,6] <- as.numeric(as.character(ggData[,6]))
ggData$sig <- as.numeric(0)
ggData$sig[(ggData$upper> 0 & ggData$lower > 0 )|(ggData$upper< 0 & ggData$lower < 0 )  ] <- 1
ggData$gov_f <- ""
ggData$gov_f[ggData$gov == 0 ] <- "Opposition"
ggData$gov_f[ggData$gov == 1 ] <- "Government"

###############################################################################
#                                Figure 4a                                     #
###############################################################################


yticks <- c(1.8,1.5,1.2)
xticks <- c(-0.4,-0.2,0)
pdf("C:/Users/rehmerjo/Dropbox/HERTIE/paper_kobo/R&R/tex/kobo_gov.pdf",width = 14, height = 6)
par(mfrow = c(1,2),
    oma = c(3,1,1,1),
    mai = c(1, 0.1, 0.1, 0.1),
    mar =  c(2, 5, 0, 1) + 0.1,
    mgp = c(3, 1, 1), 
    xpd = NA,
    bty = 'n')
plot(ggData$ame[ggData$ci == 95 & ggData$iv == "kobo" & ggData$gov_f == "Government"], yticks, 
     ylab = "", xlab = "", pch = 19, xlim = c(-0.5,0), ylim = c(1.1,1.9), 
     cex.lab = 1.6, frame = FALSE, xaxt='n', yaxt = 'n', main = "Government", cex.main = 1.75)
segments( ggData$lower[ggData$iv == "kobo" & ggData$ci == 95 & ggData$gov_f == "Government"], yticks,
          ggData$upper[ggData$iv == "kobo" & ggData$ci == 95 & ggData$gov_f == "Government"] , yticks, lty = 2   )
segments( ggData$lower[ggData$iv == "kobo" & ggData$ci == 90 & ggData$gov_f == "Government"], yticks,
          ggData$upper[ggData$iv == "kobo" & ggData$ci == 90 & ggData$gov_f == "Government"] , yticks, lty = 1   )
text(x = rep(-0.5,3), y = yticks, labels = c("PQs","PMBs init.","PMBs cosp."), cex = 1.5)
lines(x = rep(0,3), y = c(1.1,1.5,1.9), col = "red", lty = "dashed")
text(par("usr")[3], x = xticks, labels = xticks, 
     pos = 3, xpd = TRUE, cex = 1.25)
axis(side = 1, at = -0.2,labels = "Average Marginal Effects", col = "white", cex.axis = 1.75, line = 0)


xticks <- c(-60,-30,0)
plot(ggData$ame[ggData$ci == 95 & ggData$iv == "kobo" & ggData$gov_f == "Opposition"], yticks, 
     ylab = "", xlab = "", pch = 19, xlim = c(-60,1), ylim = c(1.1,1.9), 
     cex.lab = 1.6, frame = FALSE, xaxt='n', yaxt = 'n' , main = "Opposition", cex.main = 1.75)
segments( ggData$lower[ggData$iv == "kobo" & ggData$ci == 95 & ggData$gov_f == "Opposition"], yticks,
          ggData$upper[ggData$iv == "kobo" & ggData$ci == 95 & ggData$gov_f == "Opposition"] , yticks, lty = 2   )
segments( ggData$lower[ggData$iv == "kobo" & ggData$ci == 90 & ggData$gov_f == "Opposition"], yticks,
          ggData$upper[ggData$iv == "kobo" & ggData$ci == 90 & ggData$gov_f == "Opposition"] , yticks, lty = 1   )

lines(x = rep(0,3), y = c(1.1,1.5,1.9), col = "red", lty = "dashed")
text(par("usr")[3], x = xticks, labels = xticks, 
     pos = 3, xpd = TRUE, cex = 1.25)
axis(side = 1, at = -30,labels = "Average Marginal Effects", col = "white", cex.axis = 1.75, line = 0)
dev.off()


###############################################################################
#                                Figure 4b                                     #
###############################################################################


yticks <- c(1.8,1.5,1.2)
xticks <- c(-0.1,0,0.1, 0.3)
pdf("C:/Users/rehmerjo/Dropbox/HERTIE/paper_kobo/R&R/tex/kobo_district_gov.pdf",width = 14, height = 6)
par(mfrow = c(1,2),
    oma = c(3,1,1,1),
    mai = c(1, 0.1, 0.1, 0.1),
    mar =  c(2, 5, 0, 1) + 0.1,
    mgp = c(3, 1, 1), 
    xpd = NA,
    bty = 'n')
plot(ggData$ame[ggData$ci == 95 & ggData$iv == "kobo_district" & ggData$gov_f == "Government"], yticks, 
     ylab = "", xlab = "", pch = 19, xlim = c(-0.2,0.3), ylim = c(1.1,1.9), 
     cex.lab = 1.6, frame = FALSE, xaxt='n', yaxt = 'n', main = "Government", cex.main = 1.75)
segments( ggData$lower[ggData$iv == "kobo_district" & ggData$ci == 95 & ggData$gov_f == "Government"], yticks,
          ggData$upper[ggData$iv == "kobo_district" & ggData$ci == 95 & ggData$gov_f == "Government"] , yticks, lty = 2   )
segments( ggData$lower[ggData$iv == "kobo_district" & ggData$ci == 90 & ggData$gov_f == "Government"], yticks,
          ggData$upper[ggData$iv == "kobo_district" & ggData$ci == 90 & ggData$gov_f == "Government"] , yticks, lty = 1   )
text(x = rep(-0.2,3), y = yticks, labels = c("PQs","PMBs init.","PMBs cosp."), cex = 1.5)
lines(x = rep(0,3), y = c(1.1,1.5,1.9), col = "red", lty = "dashed")
text(par("usr")[3], x = xticks, labels = xticks, 
     pos = 3, xpd = TRUE, cex = 1.25)
axis(side = 1, at = 0.1,labels = "Average Marginal Effects", col = "white", cex.axis = 1.75, line = 0)


xticks <- c(-25,0,25,50)
plot(ggData$ame[ggData$ci == 95 & ggData$iv == "kobo_district" & ggData$gov_f == "Opposition"], yticks, 
     ylab = "", xlab = "", pch = 19, xlim = c(-25,50), ylim = c(1.1,1.9), 
     cex.lab = 1.6, frame = FALSE, xaxt='n', yaxt = 'n' , main = "Opposition", cex.main = 1.75)
segments( ggData$lower[ggData$iv == "kobo_district" & ggData$ci == 95 & ggData$gov_f == "Opposition"], yticks,
          ggData$upper[ggData$iv == "kobo_district" & ggData$ci == 95 & ggData$gov_f == "Opposition"] , yticks, lty = 2   )
segments( ggData$lower[ggData$iv == "kobo_district" & ggData$ci == 90 & ggData$gov_f == "Opposition"], yticks,
          ggData$upper[ggData$iv == "kobo_district" & ggData$ci == 90 & ggData$gov_f == "Opposition"] , yticks, lty = 1   )

lines(x = rep(0,3), y = c(1.1,1.5,1.9), col = "red", lty = "dashed")
text(par("usr")[3], x = xticks, labels = xticks, 
     pos = 3, xpd = TRUE, cex = 1.25)
axis(side = 1, at = 12.5,labels = "Average Marginal Effects", col = "white", cex.axis = 1.75, line = 0)
dev.off()

